In [31]:
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import pypandoc
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from sklearn import linear_model

import matplotlib.pyplot as plt
from matplotlib import rcParams
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.nonparametric.smoothers_lowess import lowess
import numpy as np
import scipy.stats as stats

from sklearn.preprocessing import PolynomialFeatures

Exercise 9

9 a) Produce scatterplot matrix for all variables

In [2]:
auto = pd.read_csv('Auto.csv')
In [3]:
fig = px.scatter_matrix(auto)
In [4]:
fig.show()

9 b) compute matrix of correlations

In [5]:
correlation_cols = set(auto.columns)
correlation_cols.remove('name')
In [6]:
corr = auto[correlation_cols].corr()
In [7]:
corr.style.background_gradient(cmap='coolwarm')
Out[7]:
cylinders mpg weight acceleration origin year displacement
cylinders 1 -0.77626 0.897017 -0.504061 -0.564972 -0.346717 0.95092
mpg -0.77626 1 -0.831739 0.422297 0.563698 0.581469 -0.804443
weight 0.897017 -0.831739 1 -0.419502 -0.581265 -0.3079 0.933104
acceleration -0.504061 0.422297 -0.419502 1 0.210084 0.282901 -0.544162
origin -0.564972 0.563698 -0.581265 0.210084 1 0.184314 -0.610664
year -0.346717 0.581469 -0.3079 0.282901 0.184314 1 -0.369804
displacement 0.95092 -0.804443 0.933104 -0.544162 -0.610664 -0.369804 1

9 c) multiple linear regression with mpg as response

In [9]:
auto['horsepower'] = pd.to_numeric(auto['horsepower'],errors='coerce')
In [10]:
auto = auto.dropna(axis=0,how='any')
In [11]:
input_cols = correlation_cols
input_cols.remove('mpg')
In [12]:
auto.shape
Out[12]:
(392, 9)
In [13]:
input_cols
Out[13]:
{'acceleration',
 'cylinders',
 'displacement',
 'horsepower',
 'origin',
 'weight',
 'year'}
In [14]:
X = auto[input_cols]
In [15]:
y = auto[['mpg']]
In [16]:
X = sm.add_constant(X)
C:\Users\clnau\AppData\Local\Continuum\anaconda3\lib\site-packages\numpy\core\fromnumeric.py:2389: FutureWarning:

Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.

In [17]:
model = sm.OLS(y,X).fit()
In [18]:
predictions = model.predict(X)
In [19]:
model.summary()
Out[19]:
OLS Regression Results
Dep. Variable: mpg R-squared: 0.821
Model: OLS Adj. R-squared: 0.818
Method: Least Squares F-statistic: 252.4
Date: Sun, 19 Jul 2020 Prob (F-statistic): 2.04e-139
Time: 21:21:45 Log-Likelihood: -1023.5
No. Observations: 392 AIC: 2063.
Df Residuals: 384 BIC: 2095.
Df Model: 7
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const -17.2184 4.644 -3.707 0.000 -26.350 -8.087
cylinders -0.4934 0.323 -1.526 0.128 -1.129 0.142
weight -0.0065 0.001 -9.929 0.000 -0.008 -0.005
horsepower -0.0170 0.014 -1.230 0.220 -0.044 0.010
acceleration 0.0806 0.099 0.815 0.415 -0.114 0.275
origin 1.4261 0.278 5.127 0.000 0.879 1.973
year 0.7508 0.051 14.729 0.000 0.651 0.851
displacement 0.0199 0.008 2.647 0.008 0.005 0.035
Omnibus: 31.906 Durbin-Watson: 1.309
Prob(Omnibus): 0.000 Jarque-Bera (JB): 53.100
Skew: 0.529 Prob(JB): 2.95e-12
Kurtosis: 4.460 Cond. No. 8.59e+04


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.59e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

9 c i) Is there a relationship between predictors and response?

Yes, there is a relationship between the predictors and the response based on the r-squared value of 0.821, which shows that 82% of the mpg response can be explained by the predictors in the model. Some predictors don't have as much statistical significance to the response though, based on their p-values.

9 c ii) Which predictors appear to have statistically significant relationship to response

Displacement, Weight, Origin and Year appear to have a statistically significant relationship to the response based on their p-values.

9 c iii) What does the coefficient for year suggest?

If all other predictors held constant, the mpg increases by 0.75 for every increase in year.

9 d) Diagnostic plots of linear regression fit

In [21]:
residuals = model.resid
fitted = model.fittedvalues
smoothed = lowess(residuals,fitted)

Residual Plot

In [22]:
fig, ax = plt.subplots()
ax.scatter(fitted, residuals, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('Residuals')
ax.set_xlabel('Fitted Values')
ax.set_title('Residuals vs. Fitted')
ax.plot([min(fitted),max(fitted)],[0,0],color = 'k',linestyle = ':', alpha = .3)
Out[22]:
[<matplotlib.lines.Line2D at 0x1ec6ad508d0>]

This plot shows a non-linear relationship between the response and the predictors

Normal Q-Q Plot

In [23]:
sorted_student_residuals = pd.Series(model.get_influence().resid_studentized_internal)
sorted_student_residuals.index = model.resid.index
sorted_student_residuals = sorted_student_residuals.sort_values(ascending = True)
In [24]:
df = pd.DataFrame(sorted_student_residuals)
df.columns = ['sorted_student_residuals']
df['theoretical_quantiles'] = stats.probplot(df['sorted_student_residuals'], dist = 'norm', fit = False)[0]
rankings = abs(df['sorted_student_residuals']).sort_values(ascending = False)

fig, ax = plt.subplots()
x = df['theoretical_quantiles']
y = df['sorted_student_residuals']
ax.scatter(x,y, edgecolor = 'k',facecolor = 'none')
ax.set_title('Normal Q-Q')
ax.set_ylabel('Standardized Residuals')
ax.set_xlabel('Theoretical Quantiles')
ax.plot([np.min([x,y]),np.max([x,y])],[np.min([x,y]),np.max([x,y])], color = 'r', ls = '--')
Out[24]:
[<matplotlib.lines.Line2D at 0x1ec700aa128>]

This plot shows that the data shows a mostly normal distribution but slightly right skewed

Scale-Location Plot

In [25]:
student_residuals = model.get_influence().resid_studentized_internal
sqrt_student_residuals = pd.Series(np.sqrt(np.abs(student_residuals)))
sqrt_student_residuals.index = model.resid.index
smoothed = lowess(sqrt_student_residuals,fitted)
In [26]:
fig, ax = plt.subplots()
ax.scatter(fitted, sqrt_student_residuals, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('$\sqrt{|Studentized \ Residuals|}$')
ax.set_xlabel('Fitted Values')
ax.set_title('Scale-Location')
ax.set_ylim(0,max(sqrt_student_residuals)+0.1)
Out[26]:
(0, 2.087937146133051)

This plot shows that the residuals are randomly spread across the range of the predictors.

Residuals vs. Leverage Plot

In [27]:
student_residuals = pd.Series(model.get_influence().resid_studentized_internal)
student_residuals.index = model.resid.index
df = pd.DataFrame(student_residuals)
df.columns = ['student_residuals']
df['leverage'] = model.get_influence().hat_matrix_diag
smoothed = lowess(df['student_residuals'],df['leverage'])
sorted_student_residuals = abs(df['student_residuals']).sort_values(ascending = False)
In [28]:
fig, ax = plt.subplots()
x = df['leverage']
y = df['student_residuals']
xpos = max(x)+max(x)*0.01  
ax.scatter(x, y, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('Studentized Residuals')
ax.set_xlabel('Leverage')
ax.set_title('Residuals vs. Leverage')
ax.set_ylim(min(y)-min(y)*0.15,max(y)+max(y)*0.15)
ax.set_xlim(-0.01,max(x)+max(x)*0.05)
plt.tight_layout()

cooksx = np.linspace(min(x), xpos, 50)
p = len(model.params)
poscooks1y = np.sqrt((p*(1-cooksx))/cooksx)
poscooks05y = np.sqrt(0.5*(p*(1-cooksx))/cooksx)
negcooks1y = -np.sqrt((p*(1-cooksx))/cooksx)
negcooks05y = -np.sqrt(0.5*(p*(1-cooksx))/cooksx)

ax.plot(cooksx,poscooks1y,label = "Cook's Distance", ls = ':', color = 'r')
ax.plot(cooksx,poscooks05y, ls = ':', color = 'r')
ax.plot(cooksx,negcooks1y, ls = ':', color = 'r')
ax.plot(cooksx,negcooks05y, ls = ':', color = 'r')
ax.plot([0,0],ax.get_ylim(), ls=":", alpha = .3, color = 'k')
ax.plot(ax.get_xlim(), [0,0], ls=":", alpha = .3, color = 'k')
ax.annotate('1.0', xy = (xpos, poscooks1y[-1]), color = 'r')
ax.annotate('0.5', xy = (xpos, poscooks05y[-1]), color = 'r')
ax.annotate('1.0', xy = (xpos, negcooks1y[-1]), color = 'r')
ax.annotate('0.5', xy = (xpos, negcooks05y[-1]), color = 'r')
ax.legend()
plt.show()

This plot shows a potential leverage at leverage > 0.175 point which is well outside the other points leverage value

Do the residual plots suggest any unusually large outliers?

The residual plots do not show many unusually large outliers, but the leverage plot does identify one observation with large leverage

9 e) Fit linear regression models with interaction effects

In [29]:
X = auto[input_cols]
In [32]:
x_interaction = PolynomialFeatures(2, interaction_only=True, include_bias=False).fit_transform(X)
In [33]:
interaction_df = pd.DataFrame(x_interaction, index=auto.index, columns = ['cylinders','displacement','horsepower','weight','acceleration','year','origin',
                                                       'cylinders:displacement','cylinders:horsepower','cylinders:weight','cylinders:acceleration',
                                                       'cylinders:year','cylinders:origin','displacement:horsepower','displacement:weight',
                                                       'displacement:acceleration','displacement:year','displacement:origin','horsepower:weight',
                                                       'horsepower:acceleration','horsepower:year','horsepower:origin','weight:acceleration',
                                                       'weight:year','weight:origin','acceleration:year','acceleration:origin','year:origin'])
In [34]:
y = auto[['mpg']]
In [35]:
interaction_df = sm.add_constant(interaction_df)
C:\Users\clnau\AppData\Local\Continuum\anaconda3\lib\site-packages\numpy\core\fromnumeric.py:2389: FutureWarning:

Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.

In [36]:
model = sm.OLS(y,interaction_df).fit()
In [37]:
predictions = model.predict(interaction_df)
In [38]:
model.summary()
Out[38]:
OLS Regression Results
Dep. Variable: mpg R-squared: 0.889
Model: OLS Adj. R-squared: 0.881
Method: Least Squares F-statistic: 104.2
Date: Sun, 19 Jul 2020 Prob (F-statistic): 4.01e-155
Time: 21:23:07 Log-Likelihood: -929.72
No. Observations: 392 AIC: 1917.
Df Residuals: 363 BIC: 2033.
Df Model: 28
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 35.4789 53.136 0.668 0.505 -69.014 139.972
cylinders 6.9886 8.248 0.847 0.397 -9.231 23.208
displacement 0.0041 0.018 0.235 0.814 -0.030 0.039
horsepower 0.5034 0.347 1.451 0.148 -0.179 1.186
weight -5.8592 2.174 -2.696 0.007 -10.134 -1.585
acceleration -20.8956 7.097 -2.944 0.003 -34.852 -6.939
year 0.6974 0.610 1.144 0.253 -0.501 1.896
origin -0.4785 0.189 -2.527 0.012 -0.851 -0.106
cylinders:displacement 0.0004 0.001 0.399 0.690 -0.001 0.002
cylinders:horsepower 0.0116 0.024 0.480 0.632 -0.036 0.059
cylinders:weight 0.2779 0.166 1.670 0.096 -0.049 0.605
cylinders:acceleration 0.4022 0.493 0.816 0.415 -0.567 1.371
cylinders:year -0.1741 0.097 -1.793 0.074 -0.365 0.017
cylinders:origin -0.0034 0.006 -0.524 0.601 -0.016 0.009
displacement:horsepower -1.968e-05 2.92e-05 -0.673 0.501 -7.72e-05 3.78e-05
displacement:weight 0.0002 0.000 1.025 0.306 -0.000 0.001
displacement:acceleration -0.0006 0.002 -0.364 0.716 -0.004 0.003
displacement:year -0.0002 0.000 -1.056 0.292 -0.001 0.000
displacement:origin 2.472e-05 1.47e-05 1.682 0.093 -4.18e-06 5.36e-05
horsepower:weight -0.0072 0.004 -1.939 0.053 -0.015 0.000
horsepower:acceleration 0.0022 0.029 0.076 0.939 -0.055 0.060
horsepower:year -0.0058 0.004 -1.482 0.139 -0.014 0.002
horsepower:origin -8.491e-05 0.000 -0.294 0.769 -0.001 0.000
weight:acceleration 0.4583 0.157 2.926 0.004 0.150 0.766
weight:year 0.0556 0.026 2.174 0.030 0.005 0.106
weight:origin -0.0035 0.003 -1.041 0.299 -0.010 0.003
acceleration:year 0.1393 0.074 1.882 0.061 -0.006 0.285
acceleration:origin 0.0240 0.019 1.232 0.219 -0.014 0.062
year:origin 0.0059 0.002 2.482 0.014 0.001 0.011
Omnibus: 47.137 Durbin-Watson: 1.664
Prob(Omnibus): 0.000 Jarque-Bera (JB): 116.060
Skew: 0.600 Prob(JB): 6.28e-26
Kurtosis: 5.380 Cond. No. 3.76e+08


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.76e+08. This might indicate that there are
strong multicollinearity or other numerical problems.
In [39]:
model.pvalues[model.pvalues < 0.05]
Out[39]:
weight                 0.007354
acceleration           0.003446
origin                 0.011921
weight:acceleration    0.003655
weight:year            0.030331
year:origin            0.013516
dtype: float64

In the interaction model above with all of the interactions included the displacement:weight, displacement:origin, and acceleration:origin interactions all appear to be significant. This shows that although origin does not appear to be significant on its own it is important in relationship to displacement and acceleration so it should be included in the model

9 f) Try a few different transformations of the variables

In [40]:
auto.columns
Out[40]:
Index(['mpg', 'cylinders', 'displacement', 'horsepower', 'weight',
       'acceleration', 'year', 'origin', 'name'],
      dtype='object')
In [43]:
auto['sqrt_displacement']=auto['displacement']**(1/2)
auto['log_displacement']=np.log(auto['displacement'])
auto['displacement_squared']=auto['displacement']**(2)
In [44]:
fig = px.scatter(auto, x='sqrt_displacement', y='mpg')
fig.show()
In [45]:
fig = px.scatter(auto, x='log_displacement', y='mpg')
fig.show()
In [46]:
fig = px.scatter(auto, x='displacement_squared', y='mpg')
fig.show()
In [47]:
fig = px.scatter(auto, x='displacement', y='mpg')
fig.show()

From the plots above for displacement, sqrt(displacment), log(displacment) and displacement^2, it appears log and square root show more linear relationship with mpg than displacement and dispalcement^2.

Exercise 10

10 a) Fit multiple regression model to predict sales using price, urban and us

In [57]:
carseats = pd.read_csv('Carseats.csv')
In [58]:
carseats['Urban'] = carseats.Urban.map(dict(Yes=1, No=0))
carseats['US'] = carseats.US.map(dict(Yes=1, No=0))
In [59]:
carseat_input_cols = ['Price', 'Urban', 'US']
In [60]:
X = carseats[carseat_input_cols]
In [61]:
y = carseats['Sales']
In [62]:
X = sm.add_constant(X)
In [63]:
model = sm.OLS(y,X).fit()
In [64]:
predictions = model.predict(X)
In [65]:
model.summary()
Out[65]:
OLS Regression Results
Dep. Variable: Sales R-squared: 0.239
Model: OLS Adj. R-squared: 0.234
Method: Least Squares F-statistic: 41.52
Date: Sun, 19 Jul 2020 Prob (F-statistic): 2.39e-23
Time: 21:23:54 Log-Likelihood: -927.66
No. Observations: 400 AIC: 1863.
Df Residuals: 396 BIC: 1879.
Df Model: 3
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 13.0435 0.651 20.036 0.000 11.764 14.323
Price -0.0545 0.005 -10.389 0.000 -0.065 -0.044
Urban -0.0219 0.272 -0.081 0.936 -0.556 0.512
US 1.2006 0.259 4.635 0.000 0.691 1.710
Omnibus: 0.676 Durbin-Watson: 1.912
Prob(Omnibus): 0.713 Jarque-Bera (JB): 0.758
Skew: 0.093 Prob(JB): 0.684
Kurtosis: 2.897 Cond. No. 628.


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

10 b) Interpretations of each coefficient in the model

When the other two predictors are held constant, when price is increased by $1, sales dicreases by 0.0545 sales. The binary (yes/no) urban predictor is not statistically significant to sales. Stores in the US sell 1200 more units than non US stores, with all other predictors held constant

10 c) Write model in equation form:

Sales = 13.0435 - 0.545 Price - 0.0219 Urban + 1.2 * US + $\epsilon$ where Urban and US are binary (Yes=1 and No=0)

10 d) For which predictors can you reject null hypothesis

We can reject null hyptohesis for Price and US, as their P-values are less than 0.05

10 e) Fit smaller model based on response to previous question

In [66]:
carseat_input_cols = ['Price', 'US']
In [67]:
X = carseats[carseat_input_cols]
In [68]:
y = carseats['Sales']
In [69]:
X = sm.add_constant(X)
In [70]:
model = sm.OLS(y,X).fit()
In [71]:
predictions = model.predict(X)
In [72]:
model.summary()
Out[72]:
OLS Regression Results
Dep. Variable: Sales R-squared: 0.239
Model: OLS Adj. R-squared: 0.235
Method: Least Squares F-statistic: 62.43
Date: Sun, 19 Jul 2020 Prob (F-statistic): 2.66e-24
Time: 21:24:03 Log-Likelihood: -927.66
No. Observations: 400 AIC: 1861.
Df Residuals: 397 BIC: 1873.
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 13.0308 0.631 20.652 0.000 11.790 14.271
Price -0.0545 0.005 -10.416 0.000 -0.065 -0.044
US 1.1996 0.258 4.641 0.000 0.692 1.708
Omnibus: 0.666 Durbin-Watson: 1.912
Prob(Omnibus): 0.717 Jarque-Bera (JB): 0.749
Skew: 0.092 Prob(JB): 0.688
Kurtosis: 2.895 Cond. No. 607.


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

10 f) How well do the models in (a) and (e) fit the data?

Both have R-squared values of 0.239, meaning they only explain 24% of the response in the prediction.

10 g) Using model from (e), obtain 95% confidence interval for the coefficients

From model summary in (e): price: [-0.065, -0.044]; US: [0.692, 1.708]

10 h) Is there evidence of outliers or high leverage observations in model from (e)?

In [73]:
residuals = model.resid
fitted = model.fittedvalues
smoothed = lowess(residuals,fitted)

Residual Plot

In [74]:
fig, ax = plt.subplots()
ax.scatter(fitted, residuals, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('Residuals')
ax.set_xlabel('Fitted Values')
ax.set_title('Residuals vs. Fitted')
ax.plot([min(fitted),max(fitted)],[0,0],color = 'k',linestyle = ':', alpha = .3)
Out[74]:
[<matplotlib.lines.Line2D at 0x1ec73a2ce48>]

Normal Q-Q Plot

In [75]:
sorted_student_residuals = pd.Series(model.get_influence().resid_studentized_internal)
sorted_student_residuals.index = model.resid.index
sorted_student_residuals = sorted_student_residuals.sort_values(ascending = True)
df = pd.DataFrame(sorted_student_residuals)
df.columns = ['sorted_student_residuals']
df['theoretical_quantiles'] = stats.probplot(df['sorted_student_residuals'], dist = 'norm', fit = False)[0]
rankings = abs(df['sorted_student_residuals']).sort_values(ascending = False)

fig, ax = plt.subplots()
x = df['theoretical_quantiles']
y = df['sorted_student_residuals']
ax.scatter(x,y, edgecolor = 'k',facecolor = 'none')
ax.set_title('Normal Q-Q')
ax.set_ylabel('Standardized Residuals')
ax.set_xlabel('Theoretical Quantiles')
ax.plot([np.min([x,y]),np.max([x,y])],[np.min([x,y]),np.max([x,y])], color = 'r', ls = '--')
Out[75]:
[<matplotlib.lines.Line2D at 0x1ec73ee6c18>]

Scale-Location Plot

In [76]:
student_residuals = model.get_influence().resid_studentized_internal
sqrt_student_residuals = pd.Series(np.sqrt(np.abs(student_residuals)))
sqrt_student_residuals.index = model.resid.index
smoothed = lowess(sqrt_student_residuals,fitted)

fig, ax = plt.subplots()
ax.scatter(fitted, sqrt_student_residuals, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('$\sqrt{|Studentized \ Residuals|}$')
ax.set_xlabel('Fitted Values')
ax.set_title('Scale-Location')
ax.set_ylim(0,max(sqrt_student_residuals)+0.1)
Out[76]:
(0, 1.792655351290514)

Residuals vs. Leverage Plot

In [77]:
student_residuals = pd.Series(model.get_influence().resid_studentized_internal)
student_residuals.index = model.resid.index
df = pd.DataFrame(student_residuals)
df.columns = ['student_residuals']
df['leverage'] = model.get_influence().hat_matrix_diag
smoothed = lowess(df['student_residuals'],df['leverage'])
sorted_student_residuals = abs(df['student_residuals']).sort_values(ascending = False)
In [78]:
fig, ax = plt.subplots()
x = df['leverage']
y = df['student_residuals']
xpos = max(x)+max(x)*0.01  
ax.scatter(x, y, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('Studentized Residuals')
ax.set_xlabel('Leverage')
ax.set_title('Residuals vs. Leverage')
ax.set_ylim(min(y)-min(y)*0.15,max(y)+max(y)*0.15)
ax.set_xlim(-0.01,max(x)+max(x)*0.05)
plt.tight_layout()

cooksx = np.linspace(min(x), xpos, 50)
p = len(model.params)
poscooks1y = np.sqrt((p*(1-cooksx))/cooksx)
poscooks05y = np.sqrt(0.5*(p*(1-cooksx))/cooksx)
negcooks1y = -np.sqrt((p*(1-cooksx))/cooksx)
negcooks05y = -np.sqrt(0.5*(p*(1-cooksx))/cooksx)

ax.plot(cooksx,poscooks1y,label = "Cook's Distance", ls = ':', color = 'r')
ax.plot(cooksx,poscooks05y, ls = ':', color = 'r')
ax.plot(cooksx,negcooks1y, ls = ':', color = 'r')
ax.plot(cooksx,negcooks05y, ls = ':', color = 'r')
ax.plot([0,0],ax.get_ylim(), ls=":", alpha = .3, color = 'k')
ax.plot(ax.get_xlim(), [0,0], ls=":", alpha = .3, color = 'k')
ax.annotate('1.0', xy = (xpos, poscooks1y[-1]), color = 'r')
ax.annotate('0.5', xy = (xpos, poscooks05y[-1]), color = 'r')
ax.annotate('1.0', xy = (xpos, negcooks1y[-1]), color = 'r')
ax.annotate('0.5', xy = (xpos, negcooks05y[-1]), color = 'r')
ax.legend()
plt.show()

From the Normal Q-Q plot there don't appear to be any outliers and from the Residuals vs. Leverage plot there don't appear to be any leverage points